Strong thermalization of a mesoscopic two-component Bose-Hubbard model 
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We study thermalization of a two-component Bose-Hubbard model by exact diagonalization. 
Initially the two components do not interact and are each at equilibrium but with different tem- 
peratures. As the on-site inter-component interaction is turned on, perfect thermalization occurs. 
Remarkably, not merely those simple "realistic" physical observables thermalize but even the density 
matrix of the whole system — the time-averaged density matrix of the system can be well approxi- 
mated by that of a canonical ensemble. A conjecture about this fact is put forward. 

PACS numbers: 05.70.Ln, 05.30.Jp, 05.30.-d 



For an isolated system, how and in which sense thermal 
equilibrium is reached from an initially non-equilibrium 
state, or even whether it can be reached or not, has long 
been a problem. A modern investigation came with the 
Fermi-Pasta-Ulam simulation as soon as the electric com- 
puter was available The surprising result was that 
the system exhibited a long-time periodic behavior with- 
out any sign of ergodicity, which was later ascribed to 
the integrability of the model in the continuum limit ^ . 
More recently, the problem revived again because of the 
possibility of using ultracold atoms to address it experi- 
mentally Q. Many different models, integrable [^-Q or 
non-integrable [6|-|9|, are investigated. For those models 
integrable, as expected, no thermalization, or at least no 
thermalization in the usual Gibbs ensemble sense is ob- 
served What is unexpected is that, even for some 
non-integrable models 6-8], thermalization does not oc- 
cur, at least at finite size. Moreover, even if thermaliza- 
tion does show up [1, 0] , it thermalizes only in the sense 
that some physical observables relax to the predicted 
values of a microcanonical/canonical ensemble — yet the 
time-averaged density matrix itself shares little feature 
with a microcanonical/canonical ensemble (an exception 
is where signature of this is observed). Therefore, 
the system thermalizes in a weak or pragmatic sense, 
since it is the few simple observables that are most ready 
to measure and thus of most concern. 

In this Letter, we investigate thermalization of the 
two-component Bose-Hubbard model. We find that this 
model, known as non-integrable, in some regime, does 
thermalize very well at a finite size. Remarkably, un- 
like previous works, it is not only some simple observ- 
ables that thermalize, but also the time-averaged density 
matrix (of the whole system) itself, which can be well 
approximated by a canonical ensemble density matrix. 
The motivation is to simulate the everyday experience 
that two objects initially at different temperatures, when 
brought in contact, equilibrate eventually. Here the two 
species of bosons act as the two objects. It is assumed 
that initially each component is at equilibrium in them- 



selves and at some finite but different temperatures and 
there is no interaction between them. Then at time t — 0, 
the inter-component interaction is switched on. The sub- 
sequent evolution is studied. 

The Hamiltonian is {h — ks — 1) Ht = Ha + Hi, + 
dtHab, where Ha,b are Hamiltonians of components a and 
h respectively. Explicitly, Ha = -^a Z]m=i(«Lam+i + 
h.c.) -I- ^^ Yl,m=i 'An'AnO'rnam and Ha -ir^ Hb a.s a ^ b. 
Here M is the number of sites and periodic boundary 
condition is assumed. The inter-component interaction 
is of the on-site type Hab = UabY.m=i (Anamb\j)m- The 
control function is defined as 0f<o — while 0t>o — 1. 
By assumption, the initial density matrix of the whole 
system is p{Q) = Pa(0) ® Pb{0), where the initial density 
matrices of the two components are {a = a,b) Pa{0) = 

'T' s"^""^" b)aO I- Here is the inverse tempera- 

ture of the component a, \f)a denotes the j-th eigenstate 
of Ha with eigenvalue i?^, and = ^Y^=\ e~''°^° is the 
partition function. The dimension of the Hilbert space 
"Hq of component a is = : with Na being 

the total atom number of component a. 

Now turn on the interaction. Denote the n-th eigen- 
state (with eigenvalue En) of the final Hamiltonian Ht>o 
as {ipn)- The density matrix at an arbitrary time later is 

formally p{t) = E^,i=i(^n|p(0)|^,)e-»(^"-^')*|^„)(V^;|, 
where D = DaDb is the dimension of the full Hilbert 
space % = 'Ha ®'Hb- At this point the time- averaged 
density matrix is defined as 



lim — 



D 



dtp{t)= (v^«ip(o)iv/)i^u)(^ii 



The operator p is of significant relevance for our purpose 
for multiple reasons. First, it is observable-free. Second, 
the time-averaged value of an arbitrary operator is given 
simply by (O) = lim^^oo 7 tr{p{t)0)dt = tr{pO). 
Third, the process of averaging over time is a process of 
relaxation in the sense that the entropy associated with p 
is definitely no less than that with the density matrix at 
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an arbitrary time, i.e., S{p) > S{p{t)) = S{p{0)). This is 
a corollary of the Klein inequality and is reasonable 
since p(0) contains all the information of p while the in- 
verse is invalid. The equality also means that p{t) will 
never be damped, and time-averaging is essential. 

We note that Ht is invariant under simultaneous trans- 
lations {am,bm) (a„+i,6„+i). Especially, Ht<o = 
Ha + Hh is invariant under the two translations in- 
dividually. This implies the conservation of quasi- 
momentum(a) (QM). The QM of component a is de- 
fined as Qa = X^fclo^ fcOfcO-fe (mod M), where aj. = 
'^/m e""^'^'^/^^o]„ is the creation operator of an atom 
in the fc-th Bloch state. Similar operators are defined 
for component b. Our strategy is then to transform to 
the QM space. We decompose the total Hilbert space H 
into AI subspaces according to the total QM q = Qa + Qb, 

= ®qi^^H^'^\ which are further decomposed accord- 
ing to the QM of the two components {qa,qb), ^^''^ = 
©^^o^i'^"-' H'f~'^''\ The Hamihonian and density 
matrix are always block-diagonal with respect to the q- 
subspaces, Ht = ®gL~o^H^''^ and pit) = ©^qV^^HO 
(l^ . In each g-subspace, generally there is no degen- 
eracy between the eigenstates {|V'n)}, therefore, the p 
in each g-subspace is simply the diagonal part of the 
initial density matrix in the {IfAn)} representation, i.e., 
(V'nlp^'-'IV'i) = Sni{ipn\p^'^H0)\'4^i)- Hcrc it is necessary to 
mention that for some quantities (e.g. (aj.afc)) studied 
below, we should have averaged over all the q-subspaces. 
However in this paper we do not bother doing so, be- 
cause the system behaves quantitatively similar in all the 
q-subspaces [3] • A single g-subspace captures the overall 
behavior well. Therefore, we shall focus on some specific 
q-subspace and take the normalization tr{p^'^^t)) = 1. 

As mentioned above, we are motivated to study the 
relaxation dynamics of the initial non-equilibrium sys- 
tem. A natural question is then how p^"?^ is like. As 
revealed by Fig. [TJ at least in some regime, it has strong 
characteristic of a canonical ensemble. In each panel, 
the occupations on the eigenstates p„ = {ipn\p^'^^\ipn) 
are plotted versus the eigenvalues It is amazing 

that most of the points cling close to a straight line 
except at the ends of the spectrum, and the straight 
line is actually the prediction of a canonical ensemble 
p'f\Pf) = e~^f^t>o /tr{e~^f^t>o) with the same energy 
as i.e., E = tr{H[%p^''y) = tr{H[^^J^{Pf)). The 
situation improves even further if the initial tempera- 
tures l/Pa^b are increased [3]. Besides Fig. [U we have 
explored the parameter space extensively. The rule of 
thumb is that when the J's and C/'s are comparable, i.e., 
when the model is far from the integrable limits, results 
similar to Fig.[T]can occur. Of course, we should mention 
that the fitting is not always so good as in Fig. [T] In the 
low temperature or unbalanced case {\Pa — Pb\ large), the 
fitting worsens. We will come back to this later. 

Figure [T] gives us an overall impression of To 
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FIG. 1: Diagonal elements of p^'^' (blue dots) and pc 
(red dots) versus the eigenvalues En- The parameters are 

{Na,Nb,M,q) = (4,4,9,1), {Ja,Jb,Ua,Ub) = (1,1,2,2), 

{I3a,^b) = (0.2,0.4), and (a) Uab = 0.5; (b) Uab = 1. The 
dimension of the g-subspace is D''' — 27, 225. The green 
lines depict the coarse-grained density of states of hI^q, just 
for reference (not corresponding to the vertical axes). The 
inverse temperature /?/ and the entropies of p'''' and pi''^ are 
shown in the inserts. 



characterize it further, we use the tools of distance and 
fidelity to study its relation to some reference density 
matrices. The three reference density matrices chosen 
are respectively the canonical ensemble one mentioned 
above, the product state (g-section actually) p^^rodiPf) — 
A/i ®ll=l e-'^i"-'"'^ (g) e-'^^-f^^'^'"' , and the initial density 
matrix p(9)(0) ^A/'affi^fl^e-^-^^'.^e-^''^^'""', where 
A/i,2 are normahzation coefficients such that tr{p^^^^j) — 
t7-(p('J) (0)) = 1. The distance and fidelity between two 
density matrices are defined as D{p,a) = \tr^J [p — a)"^ 
and F{p, a) = tr\J p^l'^a p^l"^ ^ respectively. They both 
take values in the range of [0, 1] and are closely related to 
each other by the inequality \ -F <D < Vl - F"^ . 
Two density matrices are close to each other if D and 1 — 
F are much smaller than unity. In Fig. [21 the Z?'s and i^'s 
are shown as the interaction strength Uoh is varied while 
all other parameters fixed. We see that in the full range 
of Uab investigated, p^c^ is the one closest to p''^, while 
p'^*-'(0) is the farthest [and so is p'^'^^{i) actually, because 
D(p(«)(0),p(«)) = D{p^'i\t),p(i^) and F(p(«)(0), p(«)) = 
F(p('?)(i),p(9') by the unitary invariance of D and F], 
with Pp^o(^ being the intermediate one. Moreover, the 

distance and infidelity between p^''^ and pi''^ are much 
smaller than unity throughout the range. This indicates 
that the "artificial" time-averaging is very effective in 
thermalizing an initially non-equilibrium state. The fact 
that p\P is always a better approximation of p^'^^ than 
Pprod indicates that the two subsystems equilibrate as a 
whole instead of factorizably. This is consistent with the 
fact that the inter-component interaction is a bulk type 
not a surface type. 

It has been verified in many aspects that p^"?^ can 
be well approximated by some canonical ensemble den- 
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FIG. 2: (a) Distance and (b) Fidelity between p^''^ and the 
canonical ensemble density matrix pi''^ (red *), product den- 
sity matrix Pp'^^ (blue o), and the initial density matrix 
p'-''\0) (black □). The parameters are {Na,Ni,M,q) = 
(4,4,9,1), {Ja,Jt,Ua,Ub) = (1,1,2,2), and (/3a, A) = 
(0.3,0.8). 



sity matrix p^c\Pf)- We then anticipate that the time- 
averaged values (predicted by p*-"^-*) of many physical 
quantities of interests are also well predicted by pc^"*. 
Moreover, if dephasing is effective, we shall see perfect 
relaxation phenomenon in the dynamics of the physical 
quantities. It is indeed the case. As shown in Fig. |31 
the occupations on the Bloch modes {a\ak) [and {h\bk)\ 
relax to their average values quickly, exhibiting minimal 
fluctuations [31 , and these values are very close to those 
predicted by p"c^ (but with much larger errors to those by 
the microcanonical ensemble density matrix p'^^^f. [l6|). 
That is, these quantities relax and relax to their equi- 
librium values in pf^. Note that due to the symmetric 



parameters chosen, {a\ak) — {b\hk) for 



and 



/o|^)^ all, and indeed we see in Fig. [3] that {a\ak) and 
{hlfik) merge from distinct initial values. Besides these 
single-particle quantities, also shown are the evolution 



of the two-particle quantities 



'^m ) , 



and a O &, 



with similar behavior observed. 

Here some remarks about the connection between ther- 
malization of the density matrix and that of physical ob- 
servables are in order. The point is that the former im- 
plies, but is not necessary for, the later. Two density ma- 
trices can yield the same expectation values for a few "re- 
alistic" physical quantities, yet be quite far apart in terms 
of D and F . Actually, it is common knowledge that in 
the thermodynamic limit, for a generic system, the pred- 
ications of a micro-canonical ensemble and a canonical 
one agree well, yet it is easy to persuade oneself that the 
distance and fidelity between the corresponding density 
matrices are (1 — D,F) <^ 1. The reason is formulated 
as the eigenstate thermalization hypothesis (ETH) [l7| . 
which is verified in some finite systems [1, Accord- 
ing to ETH, the expectation value of a generic few-body 
physical quantity varies little between eigenstates close 
in energy, therefore, the detailed distribution p„ does not 
matter as long as it is narrow in energy. Here it is veri- 
fied that ETH is acceptable for the variables in Fig. [3] (see 



FIG. 3: (a) Evolution of the occupations on the Bloch states, 
{a\ak) (solid lines) and (6].6fe) (dotted lines). From up to 
down, k = 0, . . . , 4. Other fc's are not shown because lines 
with k and M ~ k are very close to each other all the 
time, (b) Evolution of the on-site atom-atom correlations, 
{al^al^amam) (solid line) and (blnbl^bmbm) (dotted line). For 
each pair of lines, the markers of the same color on the right 
hand side indicate the average value or p prediction (*), pi''' 
prediction (□), and p^'^, prediction (o), respectively. The 
parameters are the same as in Fig. [2] with Uab = 1- 




FIG. 4: (a) Occupations on the Bloch states (one-particle op- 
erators) {a\ak) and (b) on-site atom-atom correlation (two- 
particle operator) {al„^a\^amam) for each eigenstate of H^'q. 
The t indicates the average energy E in Fig. O The parame- 
ters are the same as in Fig. (3] 



Fig.|3]). However, it plays a marginal role in the thermal- 
ization there. As shown in Fig. |31 the average energy E 
falls at the head of the spectra where ETH is not so good. 
Thus we see in Fig. [3] that the predictions of p|^)^ deviate 
significantly from the true values, yet the predictions of 
p'f^ agree much better with p^'>\ The situation persists 
in a wide range of parameters even if E falls in the body 
of the spectra where ETH is good. It is thus apparent 
that it is the detailed distribution, which is more accu- 
rately captured by p^c^ than p^]^, that really matters. 
However, this does not rule out the possibility that in 
the thermodynamic limit, the distribution of E falls in a 
small interval where ETH holds and thus both pi''' and 



P 



(-3) 



agree with p^'^' as for the physical observables. 
We now return to Fig. [TJ The fact that 
the occupations on the eigenstates p„ = 

where is a 

normalization factor, are well captured by the formula 
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Pn oc e '^J^'"", is too remarkable to be overlooked 
This fact is non-trivial since at the first glance there is no 
clue in the expression. So far we have not understood it 
fully but we do understand the weak fact that Pn/Pm — 1 
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FIG. 5: (a) Distances between the coarse-grained probability distributions corresponding to adjacent eigenstates, with Gaussion 
broadening. The green, red, and blue dots correspond to w = 0, 0.5, and 1, respectively, (b) Histogram of || P„ — Pn+i \\ in bins 
of length 0.02. (c)-(f) illustrate the Gaussion broadened {w = 0.5) probability distributions corresponding to four successive 
eigenstates. The distances between adjacent pairs are (0.0460, 0.0419, 0.0502). The parameters are the same as in Fig.[T}3 (but 
note that /3a, 6 are irrelevant here). 



if \En - E,n\ < 1//?/, in particular, Pn/p,i+i — 1. We 
have p„ = U dEadEbPn{Ea,Eb)e'^-^^--^^^^, 
with the probability distribution function Pn{Ea, Eh) — 
\{iMn)\^SiEa " E^ Eb - Ef,) [19]. It consists of a 
series of (5-functions with fixed positions but n-dependent 
amplitudes and is an intrinsic property of \ijjn) in terms 
of Coarse-graining P„ by replacing the i5-functions 
with some regular peaked function f{x,y) (satisfying 
J J dxdyf = 1 and / > 0), we rewrite p„ as 

P^^-^JJ dEadEbPn{Ea,Eb)e~^'''-~^^''K (1) 

Here the coarse-grained distribution Pn{Ea,Eb) = 
E\^J) \{ij\^n)\^fiEa - E^, Eb - E^ , and the constant 
c = J J dxdye'^^'^^^'^'''^ f{x,y), which is n-independent. 
The fact that p'^*-' and pi'^^ in Fig. [T] always agree well in 
the high temperature regime (3a,b < 0.4 suggests that 
there exists some / such that for most n's, Pn and 
Pn+i are close to each other in a certain sense — an in- 
trinsic property independent of Pa,b- As a measure of 
the difference between two probability distributions, we 
have the metric || P„ - ||= ^ // dEadEb\Pn - Pm\ 
By this metric, two distributions are close to each 
other if II P„ — Pm ||^ 1. We have studied the dis- 
tances between P„ and P„+i using the Gaussion function 
f{x,y) = ;;^exp[— (a;^ +y'^)/w'^], where w is the ad- 
justable width. The results are shown in Fig. [5] We see 
that although for the initial distributions {w = 0, in this 
case Pn degenerates to P„), || P„ — P„+i || centers around 
0.63, once broadening is triggered, it shrinks abruptly. 
For w = 0.5 and 1, over 84% and 93% pairs have a dis- 
tance less than 0.1 respectively. Moreover, those pairs 
having large distances fall mainly at the ends of the spec- 
trum, consistent with the fact that in Fig. [1] the fitting 
is bad at the ends (large fluctuations). In Figs. [5j:-f, 
the broadened distributions Pn{Ea, Eb) for four succes- 
sive n's are illustrated. It is apparent that they agree 
even in details. We also observe that the contour of P„ 
stretches along the direction Ea+Eb = const. This is rea- 
sonable since Hab, as a perturbation, mixes eigenstates of 



Eta + Hb with adjacent eigenenergies best. At this point, 
we can understand why the fitting in Fig. [T] is good and 
why low temperature and large difference |/3a — Pb\ are 
two adverse conditions for the fitting. The exponential 
weight function e~^^^'^~^''^'' descends fastest in the di- 
rection (Pa, Pb)- In the two adverse conditions, the weight 
function changes significantly across the region where P„ 
takes significant values and which extends primarily in 
the direction (1,-1). This non-uniformity potentially 
will spoil the closeness between P„ and Pn+i in terms of 
the metric above. 

Finally, we should mention that the essence of coarse- 
graining is to smear out the irrelevant details of the dis- 
tribution P„ retaining only the relevant overall informa- 
tion. As the size of the system grows, the number of 
(5-functions within the radius w increases exponentially 
and the coarse-graining shall be even more effective in 
reducing the distance between P„ and Pn-i-i, and there- 
fore it is legitimate to expect even better fitting. 

To conclude, it is demonstrated that the generic two- 
component Bose-Hubbard model exhibits perfect ther- 
malization. It is strong thermalization in that not merely 
the average values of a few physical variables but even 
the time-averaged density matrix itself thermalizes. We 
also note that our scenario is potentially realizable with 
cold atoms in optical lattices. Atoms in optical lat- 
tices are at finite temperatures necessarily [20[ and the 
inter-component interaction can be controlled with Fes- 
hbach resonance. Moreover, the ensemble average is of- 
fered automatically by a two-dimensional array of one- 
dimensional optical lattices. 

We are grateful to L. M. Duan for helpful suggestions. 
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